# clear working environment
rm(list = ls())

#### load the required packages ####
library(tidyverse)
library(PanelMatch)
library(countrycode)

#### load the clean data ####
df <- read_rds("data/dem-transitions-replication-data.rds") %>%
  glimpse()

#### make df into a properly formatted dataframe for PanelMatch ####

# begin coercing df into m
# df needs to be a dataframe class
m <- as.data.frame(df)
# check to see if m is a dataframe
class(m)

# country needs to be an integer to use PanelMatch
m$country <- as.integer(m$country)

# year needs to be an integer to use PanelMatch
m$year <- as.integer(m$year)

# create treatment variable
m$tr <- as.numeric(ifelse(m$v2x_regime >= 2, 1, 0))

# all other variables, except categorical variables, need to be numeric
m$mid <- as.numeric(m$mid)
m$fatal <- as.numeric(m$fatal)
m$ln_cinc <- as.numeric(m$ln_cinc)
m$ln_pecpp <- as.numeric(m$ln_pecpp)
m$terr_dispute <- as.numeric(m$terr_dispute)

# make sure all of the rows are unique
m <- distinct(m, country, year, .keep_all = TRUE)

# create Figure 4 in the main text
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "mid")

DisplayTreatment(panel.data = panel, legend.position = "bottom",
                        xlab = "Year", ylab = "Country",
                        title = "", legend.labels = c("Autocracy", "Democracy"),
                        color.of.treated = "black",
                        color.of.untreated = "darkgray") + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle=0, size = 6.5, vjust=0.5),
        legend.title = element_blank(), legend.key.size = unit(.25, 'cm')) + 
  scale_x_discrete(breaks = c(1925, 1950, 1975, 2000))

# save plot
ggsave("figures/treatment-plot.pdf", width = 5, height = 5, units = "in")

#### effect of democratic transitions on MID participation ####
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "mid")

# create matched sets
PMresults.CBPSweight.5.5.mid.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                               match.missing = TRUE, exact.match.variables = c("un_region"),
                                               covs.formula = ~ I(lag(mid, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                               qoi = "att", lead = 1:5, forbid.treatment.reversal = FALSE,
                                               use.diagonal.variance.matrix = TRUE)

# estimate the ATT
PEresults.CBPSweight.5.5.mid.att <- PanelEstimate(panel.data = panel, 
                                                  sets = PMresults.CBPSweight.5.5.mid.att,
                                                  se.method = "unconditional",
                                                  confidence.level = .90)

# create Figure C.1 for the appendix
msets <- PMresults.CBPSweight.5.5.mid.att
msets[["att"]][["160.1984"]]

example <- m %>%
  filter(country == 160 | country == 40 | country == 41 | country == 70 | country == 90 | country == 91 | country == 92 | country == 93 | country == 95 | country == 100 | country == 110 | country == 140 | country == 145 | country == 150 | country == 155 | country == 165) %>%
  filter(year >= 1979 & year <= 1989) %>%
  glimpse()

example <- PanelData(panel.data = example,
                     unit.id = "country",
                     time.id = "year",
                     treatment = "tr",
                     outcome = "mid")

DisplayTreatment(panel.data = example, legend.position = "bottom",
                 title = "", legend.labels = c("Autocracy", "Democracy"),
                 color.of.treated = "black",
                 color.of.untreated = "darkgray") + 
  theme(axis.text.x = element_text(angle=0, size = 6.5, vjust=0.5),
        legend.title = element_blank(), legend.key.size = unit(.25, 'cm'))  +
  labs(y = "COW Numeric Code",
       x = "Year",
       color = "Regime Type")

# save plot
ggsave("figures/example-mset.pdf", width = 5, height = 5, units = "in")

#### effect of democratic transitions on fatal MID participation ####
# create panel data
panel <- PanelData(panel.data = m,
                   unit.id = "country",
                   time.id = "year",
                   treatment = "tr",
                   outcome = "fatal")

# create matched sets
PMresults.CBPSweight.5.5.fatal.att <- PanelMatch(panel.data = panel, lag = 5, refinement.method = "CBPS.weight",
                                                match.missing = TRUE, exact.match.variables = c("un_region"),
                                                covs.formula = ~ I(lag(fatal, 1:5)) + I(lag(ln_cinc, 1:5)) + I(lag(ln_pecpp, 1:5)) + I(lag(terr_dispute, 1:5)),
                                                qoi = "att", lead = 1:5, forbid.treatment.reversal = FALSE,
                                                use.diagonal.variance.matrix = TRUE)

# estimate the ATT
PEresults.CBPSweight.5.5.fatal.att <- PanelEstimate(panel.data = panel, 
                                                    sets = PMresults.CBPSweight.5.5.fatal.att, 
                                                    se.method = "unconditional",
                                                    confidence.level = .90)

#### create a coefficients plot ####
years <- c("t+1", "t+2", "t+3", "t+4", "t+5")
estimates <- PEresults.CBPSweight.5.5.mid.att[["estimate"]]
st_error <- PEresults.CBPSweight.5.5.mid.att[["standard.error"]]

mid_gg <- data.frame(years, estimates, st_error) %>%
  mutate(lo_ci = estimates - (1.645*st_error)) %>%
  mutate(hi_ci = estimates + (1.645*st_error)) %>% 
  mutate(significant = as.factor(ifelse(lo_ci > 0 | hi_ci < 0, 1, 0))) %>%
  mutate(outcome = as.factor("All MIDs")) %>%
  glimpse() 

estimates <- PEresults.CBPSweight.5.5.fatal.att[["estimate"]]
st_error <- PEresults.CBPSweight.5.5.fatal.att[["standard.error"]]

fatal_gg <- data.frame(years, estimates, st_error) %>%
  mutate(lo_ci = estimates - (1.645*st_error)) %>%
  mutate(hi_ci = estimates + (1.645*st_error)) %>% 
  mutate(significant = as.factor(ifelse(lo_ci > 0 | hi_ci < 0, 1, 0))) %>%
  mutate(outcome = as.factor("Fatal MIDs")) %>%
  glimpse()

gg <- rbind(mid_gg, fatal_gg)

# create Figure D.1 for the Appendix
ggplot(data = gg, aes(x = years)) +
  geom_hline(aes(yintercept = 0), linetype = "dashed", linewidth = .25) +
  geom_pointrange(aes(y = estimates, ymin = lo_ci, ymax = hi_ci, color = outcome, shape = outcome), size = .4, linewidth = .70, position = position_dodge(width = .2)) +
  scale_color_manual(values = c("black", "darkgray")) +
  ylim(-.35, .35) +
  theme_minimal() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(color = "black", size = 10),
        axis.title.y = element_text(size = 10),
        legend.title = element_text(size = 10)) +
  labs(y = "Effect of Democratic Transition on Pr(Conflict)",
       x = "Years Since Transition",
       color = "",
       shape = "")

# save plot 
ggsave("figures/mid-att-plot.pdf", width = 5.5, height = 4, units = "in", bg = "white")



